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Numerical Simulation of Separated Flows over 
Arbitrary Airfoils and Their Resulting Wakes 

Introduction 

The aerodynamic design of a flight vehicle must carefully account for the 
drag. The estimation of the drag is greatly affected by viscous effects. For 
flows of practical interest, the Reynolds number is sufficiently large for the 
flow field to be divided into viscous and inviscid zones, e.g., the problem of 
flow past wing. Different approaches are available for solving such a prob- 
lem. Inherently, the Navier-Stokes formulations lead to an extremely stiff 
nonlinear system. Using an explicit algorithm to solve such problems results 
in the requirement of very small time-steps in order to satisfy the stability 
bounds. Therefore, many iterations and large computer times are required to 
reach the steady-state. To remove the time-step restriction, fully implicit 
methods have been investigated. The implicit methods, however, still require 
many iterations to reach the steady state and consequently, still require 
large computational costs. 

In an effort to decrease the computational costs associated with the im- 
plicit algorithms, many different procedures have been studied, in particular 
the Pulliam-Chaussee diagonal izat ion procedure^ and the Barth-Steger matrix 
reduction method.^ New dissipation models^ and spatially varying time-steps 
have dramatically increased the convergence rate. However, one problem still 
remains: long running times for general configuration. In aircraft design, 

any pertinent parameters must be accurately predicted (C L , C D , etc.). To this 

a 

end, high resolution is required in order to accurately compute the flow 
physics of shock and boundary-layer interaction, massive separation and turbu- 
lent flow structures. 

To overcome the problem of grid generation for complicated geometries and 
long tunning times, zonal approaches have become increasingly popular. By 
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zonal approach we mean partitioning of the flow field into distinct zones each 
of which is solved independently, where the length scales associated with each 
individual region are honored. There are a number of advantages for the zonal 
technique. First, the difficulty in generating three-dimensional grids for 
different types of complex configurations can be eliminated with the use of 
zonal methods. Second, zonal methods would allow different types of grid 
topologies to be implemented where appropriate in order for the grids to be 
mesh-efficient, that is, more points on the configuration, where accuracy is 
desired, and fewer points in the outer flow field. And finally, finer meshes 
can be used in those regions of rapid changes in the flow quantities, for ex- 
ample, in the regions where shocks occur, in the viscous boundary layer, or 
where vorticity is generated. The zonal concept has been successfully applied 
by the author 4 to some model problems for two-dimensional and axi symmetric 
flows. 

The present work is a generalization and improvement of an earlier work 
developed for studying separated flows using boundary layer type equations. 

The improvements include extensions to a general coordinate system and use of 
a more general zonal technique for solving the coupled equations. In order to 
be able to consider arbitrary geometries, second order accurate (in space) 
conservative differences are generated by considering the integral formulation 
of the governing equations in a general coordinate system. The general coor- 
dinate system is handled in as general a manner as possible to allow for the 
use of either analytically or numerically generated coordinate systems. 

The present work used a marching procedure for solving the PPNS (Par- 
tially Parabolized Navier-Stokes ) equations in the viscous region coupled in a 
fully implicit manner with the elliptic inviscid equation. To test the algo- 
rithm and compare to other solutions, solutions for flow over a flat plate and 
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flow past the symmetrical 12-percent-thick Joukowski airfoil (J012) at zero 
angle of attack were obtained. 


Analysis 

The basic equations which describe the motion of laminar incompressible 


flow are the Navier-Stokes equations. These equations can be written in 


stream function-vorticity (i|>-a>) form in general coordinates. 



f “ ii A 
'j 35 J an 1 


, 3 ,£ li . 1 it) 

3n 'j 3n j 35 ; 


( 1 ) 


and 

3(0 3co 3jj> 3u) _ _1_ , 3 .at 3^o 3. 3u> . 3_ .a jho _ J5 3(0 . 

3t 3n 35 " 35 3n R e *35 'j 35 " J 3n 3n j 3n " j 35 j 


where 

a 



(3) 


6 


x_x + y_y 

5 n 5 n 


(4) 


a 




(5) 


J = x„y - x y_ (6) 
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and R e is the Reynolds number. 

If the chosen coordinate system is orthogonal, then 0 is zero. 

Defining the unit vectors in the (£,n) coordinate system as (e^e^) 9 the 
velocity vector % can be expressed as 
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where u and v are the components of $ defined as. 
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The relation between the Cartesian velocity components (u_,v_) and the 
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present velocity components (u,v) is. 
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Evaluation of the Pressure 

The momentum equations can be written as 

p$ * v$ = -7p - yv2$ ( 12 ) 

where p is the pressure and p is the density. 

Multiplying the x-momentum equation by dx and the y-momentum equation by dy 
and adding the two to get a single equation for the pressure. 
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Equation (13) can be written as. 


<Hp ~ +P] + pu> di|> = p(|^ dy - dx) 


( 14 ) 


The Cartesian operators can be expresed in terms of the general coordinates as 
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Equation (13) is expressed in terms of the genetal coordinates as, 

dlp T- - -P-* ♦ n - f + »<- ! I? * f 


Or in terms of Cp (pressure coefficient). 
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Inviscid Analysis 

The inviscid solution in the present work is obtained from an incompress- 
ible stream function representation of the inviscid flow region with zero vor- 
ticity. The boundary conditions for the inviscid region are, at the inflow 
boundary, 5=^^, u=u w • On the interface i|> is known from the coupling between 
the viscous and inviscid zones, and for r\+<», u -ki . At the outer flow boun- 

oo 

dary, £=£ 0 , =0 (see reference 6 for more details). 

A second order accurate conservative difference scheme is generated for 
the stream function equation, by integrating the equation around a differential 
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element in the physical domain. The resulting algebraic system of equations 
is solved iteratively with the SLOR (Successive Line Over Relaxation) scheme. 
At each £-line in the computational domain, the finite difference equation at 
each nodal point is 

a j { *j.i + b A + c j { Vi = d j (19) 


where the index J denotes the grid position in the r\ direction and 5^ the 
change of between successive iterations, that is 





( 20 ) 


Along each £-line, a tridiagonal system of equations is solved using the 
Thomas algorithm. 


Viscous Analysis 

The flow in the viscous region is assumed to be governed by PPNS or (TL) 
equations , 
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The above equations are to be solved in the viscous region. 

Boundary conditions for the viscous equations are at the surface, 

q = 0 , \|; = = 0 (23a) 

at the interface. 
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r\ = ^IN' ^ “ ^IN (° btaine< * fr o m the coupling) , w * 0 (23b) 

at the inflow boundary, 

5 * a) = (i) (n) (23c) 

at the outflow boundary, 

5 = 5q# viscid/inviscid interaction is negligible. 

Central finite difference approximations are used everywhere except for 
the term (i|>^6a>£), which is treated as an upwind difference. That means, in 
the limit of the steady state, a second order accurate solution is obtained. 

At each S line, the correction equations based on a Newton linearization 
procedure have the general form: 
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A block (2x2) tridiagonal Thomas algorithm is used, where in the forward 
pass, the coefficients are calculated starting with the interface boundary 
conditions. The boundary condition at the surface and the coupling between 
viscous and inviscid zones are treated between viscous and inviscid zones are 
treated in the same way described in reference 6. 


Numerical Generation of Metric Coefficients 
As explained, one way of generating a second order accurate conservative 
differencing scheme includes calculating the metric coefficients and the 
Jacobian of the coordinate transformation at the center (in the transformed 
plane) of each side of the differential element shown in sketch 1. 
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Sketch 1, Differential Element 


On sides 1 and 3 the coefficients % and ^ will be calculated while % and 

J 


1 

J 


J J J 

will be calculated on sides 2 and 4. Using side 1 as an example, , 0^ and 

J-j were previously defined as, 

(25) 
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Along side 1, Ari = 0, therefore. 


Ax 1 = x 1 ? A ^1 
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Using the subscripts in sketch 1, 
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where A? 1 = A£ 3 = AC = constant, 
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These central differences are second order accurate in space* Therefore, 
can be represented to second order accuracy as 
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The coefficients 0 and J include the terms x^ and y^. x^ and are cal- 

culated in the same manner as x,_ and y, , and their second order accurate 
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finite difference representations are 
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Therefore, the second order accurate representations of and J 1 are 
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The finite difference expressions for the metric coefficients and 
Jacobian of the transformation at the center of the remaining element sides 
are derived analogously. Therefore, to generate the metric coefficients and 
Jacobians throughout the solution plane with central differences, the 
Cartesian coordinates of the algorithm solution points, the corners of the 
differential elements, and the points outside the boundaries, must be known. 


Results and Discussion 

The algorithm was tested for two cases. The first case was flow over a 
flat plate and the second one was flow past the symmetrical Joukowski airfoil 
J0 12 at zero angle of attack. The flat plate solution exactly matched the 
known Blasius solution. For the case of J0 12 airfoil, the choice of the grid 
is very critical for viscous computations over airfoils. There has been 
general agreement that C type grids are the best for handling the trailing 
edge and wakes for airfoils. One clear advantage of the present formulation 
is the complete independence on the way by which the grid is generated. That 
means the present scheme is capable of solving any type of grids. A C type 
grid for J012 airfoil used by NASA Langley, see figure 1, was implemented for 
the present work. Only half of the grid (Figure 2) was used for the actual 
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computations due to the symmetry of the problem. The grid shown in Figure 1 

is a coarse grid for high Reynolds number flows. Only two Reynolds number 

solutions (R =1000 and R =2000) were obtained for the grid shown in figure 2. 
e e 

Figure 3 shows the friction coefficient (c f ) over the surface of the airfoil 

at R =1000. The value of c in the vicinity of the leading edge of the 
® f 

Q 

airfoil is in complete agreement with the stagnation point flow solution 

(A =1.2326). The values of c near the trailing edge of the airfoil were 
On w f 

magnified on a larger scale to be able to see the solution in that region 
(Figure 4). As the R e increased up to 2000, the flow separated at 64% of the 
chord in accordance of the experimental data of reference 8, and the solution 
of the integral form of the boundary layer equations, see figure 5. The 
values of in the separated region are displayed on a larger scale, figure 
6. The streamline contours are given in figure 7. 

Conclusion 

The contribution of the current work lies in the generalization of the 
zonal technique for solving separated flows over any arbitrary airfoil. A 
numerical study of the aerodynamic characteristics for airfoils near stall is 
viable using the developed method. 
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FIGURE 1. C TYPE GRID FOR J012 AIRFOIL 



FIGURE 2. HALF OF THE GRID USED FOR ACTUAL COMPUTATIONS 



Figure 3. Shear Stress 



Figure 4. Detailed Study of Shear Stress Near T.E. 
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Figure 5. Shear Stress Distribution Over the J012 Airfoil. 
Re = 2000 




